R Notebook

Import Files using GUI or read.csv (use help)

data1 <- read.csv("FOTsubject1.csv", header=FALSE)
data2 <- read.csv("FOTsubject2.csv", header=FALSE)

# Extract columns: V1=time, V2=pressure, V3=flow
time1 <- data1$V1
pressure1 <- data1$V2
flow1 <- data1$V3

time2 <- data2$V1
pressure2 <- data2$V2
flow2 <- data2$V3

Select the proper variables to plot the input singal vs time. Fill the blanck denoted by ????

plot(time1, pressure1, type='l', col='red', 
     xlab='Time (s)', 
     ylab='Pressure (cmH2O)',  
     main='Excitation Signal: Pressure vs Time')

# Add grid for better visualization
grid()

COmpare in a plot the flow signal for both patients.

plot(time1, flow1, type='l', col='red', 
     xlab='Time (s)', 
     ylab='Flow (L/s)',  
     main='Flow Signal Comparison: Patient 1 (red) vs Patient 2 (blue)')
lines(time2, flow2, type='l', col='blue')

legend("topright", legend=c("Patient 1", "Patient 2"), 
       col=c("red", "blue"), lty=1)
grid()

Calculate the Forier coefficients input and output

Initialization

# Initialization 
ckp1 <- rep(0, times=61)  # Fourier coefficients pressure signal patient 1
ckp2 <- ckp1              # Fourier coefficients pressure signal patient 2
ckf1 <- ckp1              # Fourier coefficients flow signal patient 1
ckf2 <- ckp1              # Fourier coefficients flow signal patient 2

We use tt for the time vector

# Use time vector (assuming same for both patients)
tt <- time1

# Frequency range from -30 to 30 Hz
kf <- seq(-30, 30)
kk <- 1

# Calculate Fourier coefficients using exponential form
for (k in kf) {
  wk <- 2 * pi * k  # Angular frequency
  expk <- complex(real=cos(wk*tt), imaginary=-sin(wk*tt))
  
  # Calculate Fourier coefficients (integral approximation using mean)
  ckp1[kk] <- mean(pressure1 * expk)
  ckp2[kk] <- mean(pressure2 * expk)
  ckf1[kk] <- mean(flow1 * expk)
  ckf2[kk] <- mean(flow2 * expk)
  
  kk <- kk + 1
}

# Plot modulus of Fourier coefficients vs Frequency
plot(kf, Mod(ckp1), type='h', col='red', 
     xlab='Frequency (Hz)', 
     ylab='|Fourier Coefficient|',
     main='Fourier Coefficients Magnitude - Pressure Signal')
grid()

The respiratory system is excited with a multisinusoid signal containing three frequencies: 5 Hz, 12 Hz, and 19 Hz. These frequencies appear as peaks in both the positive and negative frequency ranges due to the exponential form of the Fourier series.

Yes, the main frequencies are identical in both input (pressure) and output (flow) signals. This occurs because the respiratory system behaves as a Linear Time-Invariant (LTI) system. In LTI systems, sinusoidal inputs produce sinusoidal outputs at the same frequencies, with only amplitude and phase modifications. The system does not generate new frequency components. This property is fundamental to the FOT technique, as it allows us to calculate the respiratory impedance independently for each excitation frequency.

Key Concepts to Understand 1. Period vs Frequency Period T₀ = 1 second → Signal repeats every second Fundamental frequency f₀ = 1/T₀ = 1 Hz We analyze harmonics: k = -30 to +30 Hz 2. Complex Fourier Coefficients Each coefficient c_k is a complex number Magnitude Mod(c_k): How strong is this frequency component? Phase Arg(c_k): What’s the time shift of this component? 3. Modulus (Magnitude) Mod(complex_number) # Returns magnitude/absolute value This tells you the “strength” or “energy” at each frequency 4. Negative Frequencies Mathematical artifact of exponential form Physically, positive and negative frequencies represent the same sinusoid We usually report only positive frequencies (5, 11, 19 Hz)

# Identify main frequencies (look for peaks)
cat("\nMain excitation frequencies (Hz):\n")

Main excitation frequencies (Hz):

main_freqs <- kf[Mod(ckp1) > max(Mod(ckp1))*0.1]
print(main_freqs)

[1] -19 -12  -5   5  12  19

Synthesize signals using main Fourier coefficients

Typically, the main frequencies are around 5, 12, 19 Hz and their negatives. We should identify these from the plot above. For this example, let’s assume main frequencies are -19, -12, -5, 5, 12, 19

ff <- c(-19, -12, -5, 5, 12, 19)

# Find indices for these frequencies
ii <- kf %in% ff

# Synthesize the pressure signal using only main coefficients
pressure_synth <- rep(0, length(tt))
for (idx in which(ii)) {
  k <- kf[idx]
  wk <- 2 * pi * k
  expk_synth <- complex(real=cos(wk*tt), imaginary=sin(wk*tt))
  pressure_synth <- pressure_synth + Re(ckp1[idx] * expk_synth)
}

# Compare empirical vs synthesized signal
plot(tt, pressure1, type='l', col='gray', 
     xlab='Time (s)', 
     ylab='Pressure (cmH2O)',
     main='Empirical (gray) vs Synthesized (red) Pressure Signal')
lines(tt, pressure_synth, col='red', lwd=2)

legend("topright", legend=c("Empirical", "Synthesized"), 
       col=c("gray", "red"), lty=1, lwd=c(1,2))
grid()

cat("\nThe synthesized signal has LESS noise because we only use the main\n")

The synthesized signal has LESS noise because we only use the main

cat("frequency components and filter out high-frequency noise.\n")

frequency components and filter out high-frequency noise.

Compute Complex Impedance and Plot Results

# Calculate respiratory impedance: Z = Pressure / Flow
# For the selected frequencies only
Zresp1 <- ckp1[ii] / ckf1[ii]
Zresp2 <- ckp2[ii] / ckf2[ii]

# Use only positive frequencies for plotting (traditional convention)
positive_freq_mask <- ff > 0
ff_positive <- ff[positive_freq_mask]
Zresp1_pos <- Zresp1[positive_freq_mask]
Zresp2_pos <- Zresp2[positive_freq_mask]

# Plot REAL part of respiratory impedance (Resistance)
plot(ff_positive, Re(Zresp1_pos), type='b', col='blue', pch=19,
     xlab='Frequency (Hz)', 
     ylab='Resistance (cmH2O·s/L)',
     main='Real Part of Respiratory Impedance',
     ylim=c(-6, 6))
points(ff_positive, Re(Zresp2_pos), type='b', col='red', pch=19)

legend("topright", legend=c("Patient 1", "Patient 2"), 
       col=c("blue", "red"), pch=19, lty=1)
grid()

# Plot IMAGINARY part of respiratory impedance (Reactance)
plot(ff_positive, Im(Zresp1_pos), type='b', col='blue', pch=19,
     xlab='Frequency (Hz)', 
     ylab='Reactance (cmH2O·s/L)',
     main='Imaginary Part of Respiratory Impedance',
     ylim=c(-6, 6))
points(ff_positive, Im(Zresp2_pos), type='b', col='red', pch=19)

abline(h=0, lty=2, col='gray')
legend("topright", legend=c("Patient 1", "Patient 2"), 
       col=c("blue", "red"), pch=19, lty=1)

grid()

From Question 5, you identified 3 main excitation frequencies, for example: 5 Hz, 12 Hz, 19 Hz Due to the exponential form symmetry, each frequency appears twice: Positive frequencies: +5, +12, +19 Hz → 3 coefficients Negative frequencies: -5, -12, -19 Hz → 3 coefficients Total: 6 coefficients

Mathematical Relationship For a real signal (like pressure):

c_k and c_{-k} are complex conjugates
c_{-k} = c_k* (the asterisk means complex conjugate)

They work together to create real sinusoids Understanding Signal Synthesis The Fourier Synthesis Formula General formula (using all coefficients):

# Initialize synthesized signal
pressure_synth <- rep(0, length(tt))

# For each main frequency
for (idx in which(ii)) {  # ii = logical vector for main frequencies
  k <- kf[idx]            # Get frequency value
  wk <- 2 * pi * k        # Angular frequency
  
  # Complex exponential: e^(j·ω·t)
  expk_synth <- complex(real=cos(wk*tt), imaginary=sin(wk*tt))
  
  # Add this frequency component to the synthesized signal
  pressure_synth <- pressure_synth + Re(ckp1[idx] * expk_synth)
}

# Full signal comparison
plot(tt, pressure1, type='l', col='gray60', lwd=1,
     xlab='Time (s)', 
     ylab='Pressure (cmH2O)',
     main='Signal Synthesis: Empirical (gray) vs Synthesized (red)')
lines(tt, pressure_synth, col='red', lwd=2)

legend("topright", 
       legend=c("Empirical (with noise)", "Synthesized (6 main frequencies)"), 
       col=c("gray60", "red"), lty=1, lwd=c(1,2))
grid()

# Zoom in on a small time window to see noise clearly
time_window <- tt >= 0 & tt <= 2  # First 2 seconds
plot(tt[time_window], pressure1[time_window], type='l', col='gray60', lwd=1,
     xlab='Time (s)', 
     ylab='Pressure (cmH2O)',
     main='Zoomed View: Noise Comparison (First 2 seconds)')
lines(tt[time_window], pressure_synth[time_window], col='red', lwd=2)

legend("topright", 
       legend=c("Empirical (noisy)", "Synthesized (smooth)"), 
       col=c("gray60", "red"), lty=1, lwd=c(1,2))
grid()

# Step 6.4: Quantify the noise
cat("\n--- STEP 4: Noise Analysis ---\n")

--- STEP 4: Noise Analysis ---

# Calculate the difference (this is the filtered noise)
noise_estimate <- pressure1 - pressure_synth

# Calculate signal power and noise power
signal_power <- mean(pressure_synth^2)
noise_power <- mean(noise_estimate^2)
snr <- 10 * log10(signal_power / noise_power)

cat(sprintf("Signal Power: %.4f\n", signal_power))

Signal Power: 1.5290

cat(sprintf("Noise Power: %.4f\n", noise_power))

Noise Power: 0.0481

cat(sprintf("Signal-to-Noise Ratio (SNR): %.2f dB\n", snr))

Signal-to-Noise Ratio (SNR): 15.02 dB

cat(sprintf("Noise RMS: %.4f cmH2O\n", sqrt(noise_power)))

Noise RMS: 0.2193 cmH2O

# Plot the extracted noise
plot(tt[time_window], noise_estimate[time_window], type='l', col='darkred',
     xlab='Time (s)', 
     ylab='Difference (cmH2O)',
     main='Extracted Noise: Empirical - Synthesized')
abline(h=0, lty=2, col='gray')

grid()

DIAGNOSTIC INTERPRETATION

cat("\n========== DIAGNOSTIC ANALYSIS ==========\n")

========== DIAGNOSTIC ANALYSIS ==========

cat("\nBased on the impedance plots:\n")

Based on the impedance plots:

cat("\nREAL PART (Resistance):\n")

REAL PART (Resistance):

cat("- Higher resistance, especially at low frequencies → airway obstruction\n")

- Higher resistance, especially at low frequencies → airway obstruction

cat("- COPD patients show increased resistance\n")

- COPD patients show increased resistance

cat("\nIMAGINARY PART (Reactance):\n")

IMAGINARY PART (Reactance):

cat("- Resonance frequency (where reactance crosses zero)\n")

- Resonance frequency (where reactance crosses zero)

cat("- Higher resonance frequency → increased tissue stiffness (COPD)\n")

- Higher resonance frequency → increased tissue stiffness (COPD)

cat("- COPD patients show resonance at higher frequencies\n")

- COPD patients show resonance at higher frequencies

# Calculate resonance frequency (zero crossing of imaginary part)
# Simple linear interpolation
find_resonance <- function(freqs, reactance) {
  for (i in 1:(length(reactance)-1)) {
    if (reactance[i] * reactance[i+1] < 0) {
      # Linear interpolation
      f_res <- freqs[i] + (freqs[i+1] - freqs[i]) * 
               (-reactance[i]) / (reactance[i+1] - reactance[i])
      return(f_res)
    }
  }
  return(NA)
}

res_freq1 <- find_resonance(ff_positive, Im(Zresp1_pos))
res_freq2 <- find_resonance(ff_positive, Im(Zresp2_pos))

cat("\nResonance Frequencies:\n")

Resonance Frequencies:

cat(sprintf("Patient 1: %.2f Hz\n", res_freq1))

Patient 1: 10.28 Hz

cat(sprintf("Patient 2: %.2f Hz\n", res_freq2))

Patient 2: 18.46 Hz

cat("\n========== CONCLUSION ==========\n")

========== CONCLUSION ==========

if (mean(Re(Zresp1_pos)) > mean(Re(Zresp2_pos))) {
  cat("Patient 1 likely has COPD (higher resistance)\n")
  cat("Patient 2 is likely healthy\n")
} else {
  cat("Patient 2 likely has COPD (higher resistance)\n")
  cat("Patient 1 is likely healthy\n")
}

Patient 2 likely has COPD (higher resistance)
Patient 1 is likely healthy